{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Image projection manipulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ogr\n",
    "import gdal\n",
    "import rasterio\n",
    "\n",
    "image_fpath = '/Users/wronk/Data/divot_detect/B04_011293_1265_XN_53S071W.tiff'\n",
    "output_fpath = '/Users/wronk/Data/divot_detect/B04_011293_1265_XN_53S071W_warp.tiff'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "metadata = {'AREA_OR_POINT': 'Area'}\n",
      "orig_proj = PROJCS[\"SimpleCylindrical Mars\",GEOGCS[\"GCS_Mars\",DATUM[\"D_Mars\",SPHEROID[\"Mars\",3396190,0]],PRIMEM[\"Reference_Meridian\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Equirectangular\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",180],PARAMETER[\"standard_parallel_1\",0],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]]\n",
      "Validation: 0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x111e9dae0> >"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "img = gdal.Open(image_fpath)\n",
    "print(f'metadata = {img.GetMetadata()}')\n",
    "\n",
    "orig_proj = img.GetProjection()\n",
    "\n",
    "print(f'orig_proj = {orig_proj}')\n",
    "\n",
    "# Two reasonable projections\n",
    "sine_proj = '''+proj=sinu +lat_0=+52d52'33.06\" +lon_0=-71d38'19.52\" +ellps=Mars +a=3396190 +b=3396190 +units=m +no_defs'''\n",
    "#eqc_proj = '''+proj=eqc +lat_ts=+52d52'33.06\" +lat_0=0 +lon_0=180 +x_0=0 +y_0=0 +a=3396190 +b=3396190 +units=m +no_defs'''\n",
    "eqc_proj = '''+proj=eqc +lat_ts=-53.21334003281581 +lat_0=0 +lon_0=180 +x_0=0 +y_0=0 +a=3396190 +b=3396190 +units=m +no_defs'''\n",
    "\n",
    "# Verify and project\n",
    "srs = ogr.osr.SpatialReference()\n",
    "srs.ImportFromProj4(eqc_proj)\n",
    "print(f'Validation: {srs.Validate()}')\n",
    "\n",
    "gdal.Warp(output_fpath, image_fpath, dstSRS=srs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Input image:\n",
      "+a=3396190 +b=3396190 +lat_0=0 +lat_ts=0 +lon_0=180 +no_defs +proj=eqc +units=m +x_0=0 +y_0=0\n",
      "Output image:\n",
      "+a=3396190 +b=3396190 +lat_0=0 +lat_ts=-53.21334003281581 +lon_0=180 +no_defs +proj=eqc +units=m +x_0=0 +y_0=0\n"
     ]
    }
   ],
   "source": [
    "src_orig = rasterio.open(image_fpath)\n",
    "print(f'Input image:\\n{src_orig.crs}')\n",
    "\n",
    "src_warp = rasterio.open(output_fpath)\n",
    "print(f'Output image:\\n{src_warp.crs}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-53.21334003281581"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "src_orig.lnglat()[1]  # Raster center"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import requests\n",
    "rgb_path = 'https://github.com/mapbox/rasterio/raw/master/tests/data/RGB.byte.tif'\n",
    "image_rgb_byte = requests.get(rgb_path, stream=True)\n",
    "image_rgba_byte = requests.get('https://github.com/mapbox/rasterio/raw/master/tests/data/RGBA.byte.tif')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "local_save_path = '/Users/wronk/Builds/divot-detect/test.tif'\n",
    "\n",
    "\n",
    "if image_rgb_byte.status_code == 200:\n",
    "    with open(local_save_path, 'wb') as f:\n",
    "        for chunk in image_rgb_byte.iter_content(1024):\n",
    "            f.write(chunk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'local_save_path' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-44-850ede0e9867>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mskimage\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0;32mwith\u001b[0m \u001b[0mrasterio\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlocal_save_path\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0msrc\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m     \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msrc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msrs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'local_save_path' is not defined"
     ]
    }
   ],
   "source": [
    "import skimage\n",
    "with rasterio.open(local_save_path) as src:\n",
    "    print(src.srs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(6423080.0, 4.9049, 0.0, -3134200.0, 0.0, -4.9049)\n",
      "\n",
      "Driver: GTiff/GeoTIFF\n",
      "\n",
      "Projection: PROJCS[\"SimpleCylindrical Mars\",GEOGCS[\"GCS_Mars\",DATUM[\"D_Mars\",SPHEROID[\"Mars\",3396190,0]],PRIMEM[\"Reference_Meridian\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Equirectangular\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",180],PARAMETER[\"standard_parallel_1\",0],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "osgeo.gdal.Band"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "img.GetDriver()\n",
    "print(img.GetGeoTransform())\n",
    "#band = img.GetRasterBand(1)\n",
    "#print(gdal.GetDataTypeName(band.DataType))\n",
    "\n",
    "print(\"\\nDriver: {}/{}\".format(img.GetDriver().ShortName,\n",
    "                             img.GetDriver().LongName))\n",
    "print(f'\\nProjection: {img.GetProjection()}')\n",
    "\n",
    "type(band)\n",
    "#img.GetProjectionRef()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "Band_GetStatistics() takes exactly 3 arguments (1 given)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-35-ba0a9f3a2ca8>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mraster\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mGetStatistics\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/.virtualenvs/mars3/lib/python3.6/site-packages/osgeo/gdal.py\u001b[0m in \u001b[0;36mGetStatistics\u001b[0;34m(self, *args)\u001b[0m\n\u001b[1;32m   2397\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mGetStatistics\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2398\u001b[0m         \u001b[0;34m\"\"\"GetStatistics(Band self, int approx_ok, int force) -> CPLErr\"\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2399\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0m_gdal\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mBand_GetStatistics\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2400\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2401\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: Band_GetStatistics() takes exactly 3 arguments (1 given)"
     ]
    }
   ],
   "source": [
    "raster.GetStatistics()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
